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ABSTRACT 

We investigate the growth or decay rate of the fundamental mode of even symmetry 
in a viscous accretion disc. This mode occurs in eccentric discs and is known to be 
potentially overstable. We determine the vertical structure of the disc and its modes, 
treating radiative energy transport in the diffusion approximation. In the limit of very 
long radial wavelength, an analytical criterion for viscous overstability is obtained, 
which involves the effective shear and bulk viscosity, the adiabatic exponent and the 
opacity law of the disc. This differs from the prediction of a two-dimensional model. 
On shorter wavelengths (a few times the disc thickness) , the criterion for overstability 
is more difficult to satisfy because of the different vertical structure of the mode. In 
a low-viscosity disc a third regime of intermediate wavelengths appears, in which the 
overstability is suppressed as the horizontal velocity perturbations develop significant 
vertical shear. We suggest that this effect determines the damping rate of eccentricity 
in protoplanetary discs, for which the long-wavelength analysis is inapplicable and 
overstability is unlikely to occur on any scale. In thinner accretion discs and in decre- 
tion discs around Be stars overstability may occur only on the longest wavelengths, 
leading to the preferential excitation of global eccentric modes. 
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. 1 INTRODUCTION 

Accretion discs may support a wide variety of wave motions. 
The basic epicyclic tendencies of orbiting particles, which re- 
sult from gravitational and inertial forces, can combine with 
acoustic and buoyancy effects to allow many different modes 
of oscillation. In some cases self-gravitational or magnetic 
forces may also be significant. 

Waves in accretion discs are usually analysed without 
regard to the shear stress that is necessarily present in or- 
der to transport angular momentum and facilitate accre- 
tion. This stress, probably of turbulent origin, is most of- 
ten modelled as deriving from an effective viscosity, which 
might be ex pected always to damp the waves. However, 

iKatd l|l97dl found that 'pulsational instability' can occur 
under certain circumstance s, by analogy wit h the 'oversta- 
bility' of stellar oscillations iEddingtonlll926|) . In particular, 
if the stress responds to wavelike perturbations in a certain 
way, a viscous overstability is possible, in which some of 
the energy being drawn by the stress from the differential 
rotation of the disc is diverted into waves of growing ampli- 
tude. This effect has long been discussed as a means of excit- 
ing disturbances in accretion discs and planetary rings (e.g . 

iKato fc FukuellQSdlBorderies. Goldreich fc TremainJlQSSl: 



Paoaloizou & StanlevI 


1986| 


iKato. Honma & Matsumotd 


19881: iPaoaloizou & Lin 


|l988 


:ISchmit & Tscharnutedll998D. 



The stresses in accretion discs and planetary rings 
are probably not accurately described as arising from an 
isotropic viscosity in the sense of the Navier-Stokes equa- 
tion. Nevertheless, the general principle that overstability 
may occur as the stress responds to wavelike perturbations 
applies to all discs, even if the details are difficult to quan- 
tify 

Waves in thin discs are often studied in a two- 
dimensional approximation in which motion in the vertical 
(2:) direction, perpendicular to the orbital plane, is absent 
and the disc remains in vertical hydrostatic balance. This 
approximation is highly convenient but can only be justified 
under exceptional circumstances and indeed eliminat es al- 
most all of the wave modes. The original analysis of iKatol 
l|l978h was in fact based on a three-dimensional disc, but 
at the time of writing of his paper the behaviour of waves 
in such discs was not fully understood. Subsequently it was 
shown that a thin disc acts as a waveguide that conducts in 
the radial direction a set of discrete modes having distinct 
vertical structure s determine d from an eigenvalue problem 
iOkazaki fc Katolll98 5: Losk3 ll98fil : f Lubow fc Pri^ll99.?l : 
iKorvcanskv fc Pringleill995i) . Some of these can be classi- 
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fied as f, p and g modes by direct analogy with stellar os- 
cillations dOgilvidlTg QSl. Generally the propagation of these 
waves is restricted to certain intervals of radius because they 
encounter turning points at Lindblad or other resonances. 

Perhaps of greatest interest among these is the f° mode, 
the fundamental mode of even symmetry about the mid- 
plane, in the classification of , 'Ogilvie (1993)- This mode is 
the closest equivalent of the inertial-acoustic wave (or 'den- 
sity wave') that appears in the analysis of a two-dimensional 
disc. Its eigenfunction involves horizontal velocity perturba- 
tions which are symmetric about the midplane and have 
no nodes in their vertical structure, but which in general 
do vary significantly with z and are accompanied by verti- 
cal velocity perturbations. This is by far the dominant mode 
excited by tidal forcing at Lindblad resonances, and its exci- 
tation, propagation and dissipation have been investigated 
bv ILubow fe Qgilvie (1998), Ogilvie & Lubow (1999) and 
iBate et al.N2002ll . The f mode is also the obv ious candidate 
for viscous overstability. However. lKatol lll978l) assumed that 
the horizontal velocities would be independent of z, which 
is not true of the P mode. A self-consistent linear calcu- 
lation of overstable modes in a three-dimensional disc has 
never been carried out, though the nonlinear behaviour of 
axisymmetric dies with full vertical structure has been sim- 
ulated by Kley et al. (1993). They report that when a is 
larger than a small value the disc exhibits global, axisym- 
metric, overstable modes confined to the outer edge of the 
computational domain. 

The question of overstability also arises in the theory 
of eccentric discs. A small eccentricity can be regarded as 
a special type of wave motion, corresponding to azimuthal 
wavenumber m = 1, propagating in a circular disc. Although 
the viscous overstability has been discussed mainly for ax- 
isymmetric waves, the local dynamics of waves of modest 
azimuthal mode number in a thin disc is equivalent to that 
of axisymmetric waves. Studies of eccentric viscous discs 
in t wo and three dimensions have indeed found overstabi l- 
ity jLvubarskii. Postnov fc ProkhorovlligQ^ lQ'gilviell200ll) . 
meaning that the eccentricity often has a negative diffusion 
coefficient and would grow spontaneously and preferentially 
on small radial scales (comparable to the disc thickness). 
This is contrary to the intuitive concept that viscosity should 
damp eccentricity. It is important to note that m = 1 waves 
propagating in a (nearly) Keplerian disk are quite excep- 
tional in that they can maintain a long radial wavelength 
over an extended radial distance. 

Understanding the growth or decay rate of eccentricity 
in a disc is of great importance. For example, in protoplan- 
etary systems, eccentricity is a joint property of the planets 
and the disc since it is exchanged through secular and reso- 
nant interactions. The behaviour of eccentricity within the 
disc therefore affects the eccentricity of all the planets as 
well. 

The decretion discs around Be stars appear to have 
global eccent ric modes which may be ex cited by the viscous 
overstability dKat dll983l:IOkazakilll99 j) . However, since the 
overstability is believed to grow preferentially on a wave- 
length comparable to the disc thickness it may be difficult 
to account for the ap pearance of coh erent global modes. 

It was shown by lOgilvid (1200 j) that the criterion for 
overstability is different in two- and three-dimensional discs, 
and that it can be suppressed either by introducing a suf- 



ficiently large bulk viscosity or by using a more sophisti- 
cated model of the turbulen t stress that invo l ves a suffi- 
ciently large relaxation time. iLatter & Qgilvie' f2006') have 
also shown that the viscous overstability is suppressed in di- 
lute planetary rings when the non-Newtonian nature of the 
stress is taken into account. The analysis of eccentric discs by 
:Ogi lvij i200 J) assumes that the eccentricity is independent 
of z, which is correct if sufficiently large stresses are present 
to couple different strata in the disc. In a low-viscosity disc, 
however, the eccentricity could vary with z and might be 
better described by the structure of the P mode. 

We are led, therefore, to investigate the behaviour of 
wave modes in a viscous disc. In Section|5|we analyse a sim- 
ple two-dimensional model by way of introduction, while the 
full three-dimensional model is presented in Section |3] We 
describe the methods of numerical solution in Section |1| and 
the results follow in Section|^ We conclude with a summary 
and discussion in Section |n] 



2 TWO-DIMENSIONAL MODEL 

In order to describe disturbances with wavelengths that are 
short compared to the r adial coordinate r we use the m odel 
of the shearing sheet jGoldreich fc Lvnden-Be^illl965^ . A 
point in the disc at reference radius ro and orbiting with 
angular velocity f2o is used as the origin of a rotating 
Cartesian coordinate system in which x, y and z corre- 
spond to the radial, azimuthal and vertical directions. The 
local shear rate in the disc is given by Gort's first con- 
stant ^0 = — (7'o/2)(df2/dr)o. Having selected the reference 
point, we subsequently omit the subscript zero. We assume 
throughout that Q, > A> Q\m & Keplerian disc, A = 

In this section, by way of introduction and for the sake 
of comparison with subsequent results, we consider a sim- 
ple two-dimensional model without self-gravity, heating or 
cooling. We work with the equation of motion, 

E(9tUi -I- UjdjUi + 2Q,€ajUj) = — E9i<E> - diP + djTij, (1) 

and the equation of mass conservation, 

dt^ + d,(T.u^) ^ 0, (2) 

where E, P, and Tij are the vertically integrated density, 
pressure and viscous stress, Ui is the fluid velocity, and "1> = 
~2QAx^ is the tidal potential. The viscous stress is given 
by 

Tij = vT.{diUj + djUi) + {uh - ^u)T,{dt:Uk)5ij, (3) 

where i> and i^b are the vertically averaged kinematic viscos- 
ity and bulk viscosity. The basic state of the disc is given by 
Uy — —2Ax, E = constant, and Try = ~2AvT,. We assume 
that P, V, and v\, are specified functions of E, and define 
the sound speed through v1 = dP/dE. 

We introduce perturbations of the form Re[E' exp(st -|- 
ifca;)], etc., where s is the (complex) growth rate and k is 
the (real) radial wavenumber. These perturbations are in- 
dependent of y and correspond to axisymmetric waves in 
cylindrical geometry, but also accurately describe the local 
properties of non-axisymmetric waves of modest azimuthal 
wavenumber in a thin disc. The linearized equations read 

E(sk:, - 20^;) = -\kvlY! + ifcT4, (4) 
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T, [su'y + 2{Q ~ A)u'^] = 
sE' + Eifcu^ = 0, 



= -24^S' 



(5) 
(6) 
(7) 

(8) 



Solution of this algebraic system yields the dispersion rela- 
tion 

+ (j^b + ^i>)k^s^ + [fir + ^B^'^ + ^'(^'b + s 

+4nA^^e + v!k'uk^ = 0, (9) 
dzj 

where 57^ — 457(f2 — yl) is the square of the radial epicyclic 
frequency, and is assumed to be positive. (The symbol k, 
conventional for the epicyclic frequency, is used below to 
denote opacity.) In a Keplerian disc, Qr ~ fi- 

In the special case of an inviscid disc we obtain either a 
zero-frequency mode with s = and u'^ = 0, corresponding 
to a steady 'geostrophic flow' in which the Coriolis force is 
balanced by a pressure gradient, or an undamped inertial- 
acoustic wave with frequency uj given by the familiar disper- 
sion relation ui^ = — — + v'^k^ . 

Consider first the long- wavelength limit, fc — > 0, of the 
dispersion relation. One root behaves according to 



k^ + 0{k*) 



(10) 



0.2 dE 

which yields exponential growth (viscous instability) when 
d(i7E) 



dE 



<0, 



as found originally by iLightman &: EardlevI lll97' 

other two roots have the behaviour 

s = ±ia- ( 1 "'^'^ 



2n^ 



1 

+ 2 



A:VlA d(i^E) 



^2 dE 



0(fc* 



(11) 
The 



(12) 



This yields exponentially growing oscillations (viscous over- 
stability) when 



dE 



'^''UOA- 



(13) 



Now consider shorter wavelengths. If viscous instability 
occurs in the limit — * 0, we reach marginal stability with 
s = when 



~4D.A 



d(£E^ 
dE ■ 



(14) 



Alternatively, if viscous overstability occurs in the limit k 
0, we reach marginal stability with s = ~iuj ^ Q when 

v{v^ + ^v){vh + lv)k'^ + {vh + ^v)vlk'^ 



= AQ,A 



3 

d(£E) 
dE 



(t'b - 



(15) 



In each case the right-hand side of the equation is positive 
while the left-hand side is a monotonically increasing func- 
tion of k^ , vanishing at k^ — 0, so we find a unique positive 
value of k^ for marginal stability. The instability or oversta- 
bility is quenched for all shorter wavelengths. 

For example, a Keplerian disc with constant kinematic 
shear viscosity v <^ v1 /Q, and no bulk viscosity exhibits 



the viscous overstability on wavelengths >9i's/n, while the 
maximum growth rate occurs for a wavelength « ISiis/fi. 
These wavelengths are several times the disc thickness, 
partly because the case of constant u is close to marginal 
stability. We presume that the reason that the overstabil- 
ity is rarely or never observed in two-dimensional numerical 
simulations that include an explicit shear viscosity is either 
that the domain is too small or the numerical method has 
an effective bulk viscosity. (Shearing-box simulations of tur- 
bulent discs are invariably too small in the radial direction 
to detect any overstability.) 

The shearing- sheet approximation is valid only for 
wavenumbers satisfying \kr\ 3> 1. For wavelengths compa- 
rable to r, not only must the effects of cylindrical geometry 
be considered but also the global structure of the disc and 
the radial boundary conditions. For our purposes, the shear- 
ing sheet correctly describes the local properties of waves in 
discs. 



3 THREE-DIMENSIONAL MODEL 

Our full treatment is based on a three-dimensional shearing 
sheet. We omit self- gravity but include viscous heating and 
radiative cooling. 



3.1 Basic equations 

We work with the equation of motion, 

p{dtUi + UjdjUi -f 2Q.€is,jUj) = —pdi^ - dip + djTij, (16) 

the equation of mass conservation, 

dtp + d,{pu,) =0, (17) 
and the thermal energy equation, 

(7 - ^y^{dtp + Uidip + 'ypdiUi) = TijdiUj - diFi, (18) 
where 

Tij = p{diUj + djUi) + (^b - \p){dkUk)5ij (19) 

is the viscous stress [p and ^b being the dynamic shear and 
bulk viscosities), $ = —2Q,Ax^ + ^Q^z^ is the tidal potential 
(fi^ being the vertical epicyclic frequency, equal to f2 in a 
Keplerian disc), 

n3 



Sup 



(20) 



is the radiative energy fiux density, and the other symbols 
have their usual meanings. We assume an ideal gas with 
equation of state p — RpT and with 7 = constant. The 
properties p, p^ and k are regarded as functions of p and T. 

3.2 Equilibrium disc 

The basic state of the disc is given by Uy = —2 Ax, p — p{z), 
T = T{z), Txy = -2Ap, and = F^{z). For hydrostatic 
equilibrium. 



dzP = pgz = -pfi^z, 

and for thermal equilibrium, 

d^F, = -iA^p. 



(21) 



(22) 



4 H. N. Latter & G. I. Ogilvie 



These conditions are solved together with the constitutive 
relations to determine the equilibrium model. We adopt the 
'zero boundary conditions' whereby p = T = at the sur- 
faces z = iiH of the disc. These are appropriate to describe 
an optically thick disc and avoid the complications of match- 
ing to an atmospheric model. The surface density is 



E = 



pAz, 



(23) 



and we define the vertically averaged kinematic viscosity u 
(and analogously h'h) via 

f fidz. (24) 

J-H 

Under convenient assumptions the problem can b e re- 
duced to a standard dimensionless form ^IOgilviell2001^ . Let 
the viscosity be given hy fi = ap/Q and fj,h = Uhp/Q, 
where a and ctb are dimensionless constants, and let k = 
CkP^T^ , where C^, X and Y are constants. In particu- 
lar, for Thomson (electron scattering) opacity, X = Y = 0. 
The method of solution for the vertical structure is de- 
scribed in Appendix ^ where it is shown that vY, oc 

Y^O-O+ZX -2Y) / (G+X -2Y) 



3.3 Linearized equations 

As in Section |21 we introduce axisymmetric perturbations 
of the form 'Re[p' [z) exp(sf -I- ifca;)], etc., and obtain the lin- 
earized equations 

p{su', - 2Q.u'y) = -ikp' + ifcT^, -f- 9,TL, (25) 

P [Su'y + - A)U'^] = ikT^y + d,Ty„ {2^) 

psu'^ = p'g^ - dzP + ifcT^z + dzT!.^, (27) 

sp + u'Ap + pA = 0, (28) 

sp + u'^dzp + 7pA = TV, (29) 
with 

A — ifcit^ -f dzu^, (30) 

AT = (7 - 1) (-2Ar^j, + T^y iku'y - ikF^ ~ , (31) 

T!,^ = 2piku^ + {pb ~ Ip)A, (32) 

Txy ~ —2Ap' + piku'y, (33) 

TL = p{iku'-, + dzu'^), (34) 

Tyz = pdzu'y, (35) 

ri, = 2pdzu', + (Mb - Ip)A, (36) 



F' 



Sup 



-ifcT', 



l6aT 
3kp 

+ 



T p p 



dzT' 



<l + X)^ + {3-Y)^ 



dzT 



(38) 
(39) 



ctp 



(40) 



Note that M represents the effects of non-adiabatic heating 
and cooling. 

In general this eighth-order system of ordinary differen- 
tial equations must be solved numerically as an eigenvalue 
problem for the growth rate s, with the wavenumber A; as a 
real parameter. The boundary conditions are that the solu- 
tion should be regular at the surfaces z = ±-ff, which are 
singular points of the differential equations. The numerically 
determined solutions s(fc) define the various branches of the 
dispersion relation of the disc. 

If equations 12511 - 1291 1 are solved with the stress per- 
turbations T'ij and the non-adiabatic term M set to zero, 
we have a self-adjoint eigenvalue problem of the kind solved 
bv lKorvcanskv &: Pringla l|l995h . Indeed the structure of an 
optically thick disc with zero boundary conditions is very 
similar to that of a polytropic disc. We refer to this as the 
'inviscid problem'. 



3.4 Integral relation 

Using the preceding set of linearized equations and carry- 
ing out various integrations by parts, we obtain the integral 
relation 



J-H 



Az 



pVll\u^f + pN'^\u'^\^ H \u'^dzP + 'yp/^f 

-H y 7P 

-A'M +2sp {\iku'^ - iA|^ -I- \\iku'^ + dzu'^\^ 

-f|9.<-iA|V|iA|2) +SMb!A|^ 



X (pfc^|Uy|^ + p\dzu'y'^ -\- 2Ap' iku'y) > dz 



where 



Qz \ dz\np - -dz Inp 



(41) 



(42) 



is the square of the Brunt- Vaisala frequency. This relation 
is closely related to equation (13) in Kley et al. (1993). In 
the absence of viscous and non-adiabatic effects, it yields a 
real integral expression for the squared frequency u)'^ = —s^ 
of the mode, which has the usual variational property asso- 
ciated with self-adjoint eigenvalue problems. More generally, 
the imaginary part of this equation yields (with s = Si +isi) 



rH 

(37) sr / p(hL|' + hLn 

J-H 



Az 



2si 



AW+s* 



2^p' iku'y 



n-A^ 

p {\iku'^ — ^Af + ^\iku'^ + dzu'^f 
+ \dzu'z-^Af + \^Af) + y^\A\ 



dz 



+ 



n-A 



^ 5P {k^\u'y\'^ + \dzU'yf) 



dz. 



(43) 
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Within the first integral on the right-hand side, the first term 
quantifies non-adiabatic effects, while the second quantifies 
those of the underlying stress variation. The second integral 
on the right-hand side represents viscous dissipation. We 
rewrite the equation as 

= N + S + D, (44) 

where N, S, and D represent the three effects. 

This relation is similar to equation (3.7) of lKatd (Il978l) 
but is considerably more general because it does not assume 
that the mode is only slightly perturbed from a solution of 
the inviscid problem, ft shows that there are two effects that 
can promote growth of the wave. One is the non-adiabatic 
term J\f, if it is appropriately correlated with the velocity di- 
vergence A. (This effect is similar to the mechanism of ther- 
mal overstability in stars.) The other is the viscosity pertur- 
bation fj.' (which derives from the stress perturbation T^y), 
if it is appropriately correlated with the azimuthal velocity 
perturbation u'y. (This is the mechanism of viscous over- 
stability.) These two potentially positive contributions must 
compete with a number of negative definite terms which cor- 
respond to the damping effect of the unperturbed viscosity 
acting on the shear that is present in the mode. It is evident 
from this competition that modes with the least complicated 
vertical structure, such as the P mode, are the most likely 
to be overstable. 

3.5 Long- wavelength behaviour 

In Appendix |n| we present an asymptotic analysis of waves 
in the long- wavelength limit, fe — > 0, in the presence of vis- 
cosity. In this limit most of the modes have a damping rate 
— Re(s) comparable to a^l. Of those that do not, one behaves 
according to equation llOI I and therefore yields the same 
criterion for viscous instability as in the two-dimensional 
theory. It does so because, to a first approximation, the 
wave does not disturb the hydrostatic and thermal balances 
present in the unperturbed disc. 

The remaining long-wavelength solutions have a disper- 
sion relation of the form 

s ^ ±iQr + S2k'^ + 0{k*), (45) 

where the constant S2 is determined by an algebraic prob- 
lem derived in Appendix^ These waves yield a criterion 
for viscous overstability, Re(s2) > 0, that differs from the 
prediction of the two-dimensional analysis. The differences 
arise because the wave has an (epicyclic) period comparable 
to the time-scales for establishing the hydrostatic and ther- 
mal balances, and the vertical structure of the disc therefore 
requires a dynamical treatment. Vertical motion, in the form 
of a 'breathing motion' with u'^ oc z, plays a significant role 
in the dynamics of this mode. The overstability criterion for 
a Keplerian disc derived by t his method ag rees exactly with 
that found by the method of lOgilvi3 (|20^), who also noted 
the importance of vertical motion and non-equilibrium. 

As an example, for a Keplerian disc with constant 
(e.g. Thomson) opacity and a, Qb ^ 1, overstability occurs 
for 

Qb/a < (-337^ + 1387 - 97)/12. (46) 

In contrast, the two-dimensional analysis (using the result 
that uT, Qc T,^^^ for the case of constant opacity) predicts 



overstability for ah/ a < 8/3. When 7 = 5/3, the three- 
dimensional disc is more susceptible to overstability, as then 
the right-hand side equals 31/9 « 3.44. 

3.6 Modelling the turbulent stress 

As the ordinary molecular viscosity is too small to be ap- 
preciable, the viscous stress is presumably of turbulent ori- 
gin. The simplest 'closure' model, and the one we have 
adopted here, is the alpha prescription, which assumes a 
Newtonian stress-strain relation and an isotropic effective 
viscosity. However, both assumptions are likely to be poor 
approximations, because the turbulent stress (probably of 
magnetohydrodynamic origin) should have a non-zero relax- 
ation time comparable to the orbital period, and the effec- 
tive viscosity may be considerably anisotropic. In p rinciple 
we could a pply a more sophisticated closure (e.g. lOgilvid 
I2OOII lioOSi) . but doing so may raise more questions, and in 
fact obscure essential physics that is model-independent. It 
hence pays to examine the simpler alpha model disc first in 
detail. In this context the bulk viscosity is deployed mainly 
as a surrogate for whatever stabilizing infiuences (e.g. stress 
relaxation) may act against overstability in reality. 



4 NUMERICAL ANALYSIS 

We seek growing modes on general radial length-scales in a 
vertically stratified gaseous disc with small viscosity. This 
task requires the solution of two boundary-value problems. 
The first determines the vertical structure of the disc equi- 
librium (Appendix ^ equations lAGHAlOl , while the sec- 
ond determines the vertical structure of the linear modes 
("Section l3.3l equations l25H4Ut . For simplicity and clarity we 
adopt a constant (e.g. Thomson) opacity model and assume 
the disc is Keplerian; thus X = Y — 0, Q^/^ = 1 and 
A/Q = 3/4. Furthermore, the adiabatic exponent, 7, is set 
equal to either 7/5 or 5/3. Hence the main parameters to 
vary are a and Qb (the shear and bulk viscosities) and k (the 
radial wavenumber). The regimes in which we are interested 
require these to be small. 

Two numerical methods were employed: the 'shooting 
method' and Chebyshev collocation. The former method was 
used to solve the nonlinear equilibrium equations. The lat- 
ter method represents the ODEs of the linear system on 
a Gauss-Lobatto grid and requires solution of a general- 
ized algebraic eigenvalue problem. As the equations possess 
a singularity at the upper boundary, both methods were 
supplemented at this point by an asymptotic analysis in 
order to determine the correct limiting behaviour. A judi- 
cious choice of dependent variables was also necessary in 
both cases. The midplane boundary conditions follow from 
symmetry considerations. For example, even modes require 
= Au^jAz = du'y/dz = dp' /dz = at z = 0. 

In addition, we set the 'shooting method' upon the lin- 
earised perturbation equations, so as to check the results of 
the Chebyshev collocation method. For modest values of a 
and k the agreement between them was good. However, for 
k and a in the parameter regime of greatest interest, the 
'shooting method' either failed to converge to the overstable 
mode or settled on an eigenpair numerically 'contaminated' 
by a nearby mode. This we blame on the stiffness of the 
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Figure 1. Non-dimensionalized pressure, p», vertical radia- 
tive flux, _Ft, and temperature, T*, as functions of non- 
dimensionalized height, 2*, for an equilibrium disc with constant 
opacity (X = Y = 0). 



equations in this regime. The collocation method did not en- 
counter this difficulty. It, however, could produce eigenpairs 
which were spurious or unconverged, and these required win- 
nowing out. Also, when examining very long waves for low 
viscosities, Chebyshev collocation suffered from numerical 
error, which often compromised accurate determination of 
the growth rate, these being typically 0{ak^). That said, 
on scales for which interesting behaviour occurred, this was 
not an issue: with twenty Gauss-Lobatto nodes the over- 
stable growth rate had converged to ten significant figures. 
We present in the following section: 

(i) the vertical structure of the equilibrium (Fig. 1); 

(ii) the behaviour of the overstable mode's growth rate as 
a function of k for different a and 7, when Qb = 0, alongside 
the components of the integral relation 1441 in each case 
(Figs 2-4); 

(iii) eigenfunction structures for two illustrative parame- 
ter regimes (Figs 5 and 6); and 

(iv) the critical ratios of Qb and a for the onset of over- 
stability (Figs 7 and 8) . 



5 RESULTS 
5.1 Equilibrium 

The characteristic structure of the equilibrium state, which 
we plot in Fig. 1, differs little from a polytropic model, ex- 
cept for the addition of the vertical heat flux, which is nec- 
essarily zero at the midplane and approaches a constant on 
the boundary. This constant and the semithickness of the 
disc serve as the 'eigenvalues' of the boundary value prob- 
lem and take the approximate values 0.34 and 2.4 for the 
case considered [X = Y = 0). 



5.2 Growth rates and eigenfunctions 

We present the solution of the non-dimensionalized versions 
of equations 1251 - 141)11 . in which s = s*f2, k = k-j^/Un, 
u[ = u'i^flUH, while p', p' , T' , and F- take the same units 
as their equilibrium counterparts (see Appendix A). Since 
H ~ 2AUh, a wavenumber fe = 1 corresponds to a wave- 
length approximately 1.3 times the full disc thickness 2H. 
The growth rate is expressed in terms of the orbital fre- 
quency, the inverse of the dynamical time-scale. The stars 
will be dropped for the rest of this section. 

Amongst the various even modes exhibited by the disc, 
our numerical method found that only one has the potential 
to grow: as expected, the viscous analogue of the inviscid 
f mode, in the interval < fe < 1. No mode of odd symme- 
try grows. 

In Fig. 2a we plot the real part of the growth rate of 
the overstable mode, scaled by fc^, for modest shear vis- 
cosity and no bulk viscosity. As fc 0, this quantity ap- 
proaches a constant value approximately equal to 0.0459, 
the same value predicted by the long-wavelength theory (cf. 
equation IB451 . On short enough scales, fc ~ 1, the mode 
ceases to grow. Qualitatively this behaviour mirrors that of 
the two-dimensional model. 

In Fig. 2b we plot the components of the integral re- 
lation 1441 . the competition of which establishes growth or 
decay on the various scales. Immediately we note that the 
non-adiabatic effects, embodied in the term N, are an order 
smaller than those of viscous dissipation, D, and stress vari- 
ation, S. On long horizontal scales N contributes to oversta- 
bility, presumably through a variant of the e-mechanism in 
stars, resulting here from the increase of viscous dissipation 
in the compressed phase of the oscillation. But on sufficiently 
small horizontal scales, thermal diffusion dominates N and 
acts to smooth out disturbances (N< 0). Overstability is due 
almost entirely to the stress variation, which competes with 
viscous damping. This competition is most likely to succeed 
on the longest wavelengths, where the mode is dominated 
by horizontal motion and can avoid vertical shear. In this 
regime, weak viscous forces are sufficient to allow the mode 
to deviate from the natural form of the inviscid mode. On 
shorter wavelengths, however, the horizontal motion neces- 
sarily develops vertical shear and is accompanied by vertical 
motion, making the stress variation less efficient in driving 
overstability. The increasing horizontal shear also acts to 
suppress the overstability. 

When a is an order smaller or less, the prevalence of 
overstability on the various scales becomes more compli- 
cated, as Figs 3 and 4 reveal. Three characteristic length- 
scales emerge, which may be divided into: 

(i) the long- wavelength limit (described in Section 3.5): 
< fe< Q < 1; 

(ii) an intermediate regime of anomalous viscous dissipa- 
tion, and, consequently, stability: < a ^ fe ^ 1; and 

(iii) a regime on shorter scales, 0<Q<Cfc~l, in which 
overstability may reappear. 

In the long-wavelength regime, the overstable mode 
minimizes the amount of vertical shear in the horizontal mo- 
tion (and hence D) by maintaining u'^ and u'y effectively con- 
stant with respect to z (Fig. 5b). But as fe increases, the vari- 
ation of it^ and u'y with z increases until, by fc > a^''^, they 
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Figure 2. (a) Real part of the growth rate of the overstable mode 
divided by fe^, as a function of k, for a = 0.1, a^, = 0, X = Y = 
and 7 = 7/5. The curve asymptotes to a value 0.0459 as fc — > 0. 
The growth rate s is in units of H, and the wavenumber k is in 
units of as described at the beginning of Section 5.2. (b) 

Components of the integral relation 1441 . S refers to the stress 
variation contribution, N to the non-adiabatic contribution, and 
D to viscous dissipation. 





Figure 3. As for Fig. 2 but with a = 0.001, a^, = 0, X = Y = 
and 7 = 7/5. Note the interval of intermediate length-scales which 
are stable. 



resemble the eigenfunctions of the inviscid problem (Fig. 6). 
This complication in the mode's vertical structure leads to 
greater viscous dissipation and, subsequently, decay. How- 
ever, once the mode has assumed the limiting form of the 
inviscid eigenfunction, the amount of shear 'saturates'. Con- 
sequently, the viscous dissipation plateaus, and on shorter 
scales overstability may appear again. The viscous stress in 
this third regime may be considered a small perturbatio n 
to the inviscid equations, along the lines of iKatd il97d) . 
the principal difference being that Kato assumed that the 



Figure 4. As for Fig. 2 but with a = 0.007, ay^ = 0, X = Y = 
and 7 = 5/3. Note that for larger 7, all intermediate wavelengths 
are stable. 



perturbation was with respect to a mode with no vertical 
shear. 

We find that the imaginary part of the growth rate is 
unaffected by the increase in viscous dissipation when a 
is small. The long-wavelength theory correctly predicts the 
0{k^) correction up until intermediate scales. This result 
implies that the group velocity of the f mode, and in par- 
ticular that determining the propagation of eccentricity in 
discs, is described accurately by the three-dimensional long- 
wavelength theory. 

Odd modes require that the horizontal components of 
the motion are zero at the midplane. This means that if 
such modes are to minimize their vertical shear, and hence 
avoid the destructive effects of viscous dissipation, u'^ and 
u'y must remain near zero throughout the disc. But doing so 
will stall the mechanism driving overstability, as it requires a 
healthy correlation between Uy and /i'. This partly explains 
why only even modes exhibit overstability. 

For fe 1, a larger 7 may quench overstability (Fig. 4), 
whereas with long wavelengths a larger 7 generally stimu- 
lates overstability, as equation 14611 shows. The latter be- 
haviour originates in 7's enhancement of the pressure per- 
turbation in the S term, which, according to the long wave- 
length theory, proceeds from its dependence on W (cf. equa- 
tion IB37II . For low viscosities, the adiabatic index con- 
trols the magnitude of the vertical breathing motion via 
r^^/'^l = (7 ~ l)/7- Thus a larger 7 produces a larger W , 
larger vertical pressure advection, and consequently a larger 
S. This effect appears more important than enhanced vis- 
cous dissipation arising from the increase in vertical motion. 

When fc ~ 1, the mode's response to variations in 7 may 
be broken into two components. Firstly, an increase in 7 in- 
creases the scale upon which N becomes negative: thermal 
diffusion begins to stabilize the mode for smaller k. And, 
in fact, when 7 is sufficiently large, the non-adiabatic term 
is always negative. This behaviour springs from the sensi- 
tivity to 7 of the radiative damping terms in J\f (cf. equa- 
tiouEB. This may arise from their dependence on T', which 
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Figure 5. The u'^ and u'^ structure of the overstable mode for a = 
0.1, k = 0.01, 7 = 7/5 and = X = Y = 0. The eigenfunction 
has been normalized so that u'^ = 1 ai the boundary. As predicted 
by the long-wavelength theory, u'^ is linear in z and is 0(fc), while 
the u'^ eigenfunction is composed of an 0(1) constant component 
and an O(fc^) varying component. 

is larger for larger 7 (while the other non-adiabatic contri- 
butions, and A, do not vary significantly). Presumably, if 
the opacity is non-constant {X, Y non-zero) a version of the 
K-mechanism of overstability in stars may arise, though we 
do not explore that here. 

Secondly, and more importantly, an increase in 7 may 
contribute to an increase in the amount of shear in the 
horizontal velocities, by analogy with the inviscid problem. 
This arises via 7's influence on the vertical stratification 
of the disc, and the enh anced buoyancy forces which ensue. 
iLubow fc Pringi3 ilQQSTl established analytically that the re- 
lated 'two-dimensional mode' in an inviscid isothermal disc 
possesses u'^ and u'y which are oc exp(A'^^/n2) = exp[(7 — 
l){z/ H)^ /'y]. More generally, the amount of shear in the 
horizontal motions of an even invisci d p or f mode increases 
with N'^ and z (see equation 13 in iKorvcanskv fc Pringlj 
[199^. 



5.3 Stability curves 

Now we let ctb vary. Since Ob has a purely stabilizing effect 
on the mode, a unique (but possibly negative) critical value 
(Qb/a)c for marginal overstability can be identified for given 
values of a, 7 and fc. The smaller the value of (ab/Q:)c, the 
more likely that the overstability is suppressed by stabilizing 
infiuences such as turbulent stress relaxation. The resulting 
curves of marginal stability we plot in Figs 7 and 8, for 
Q = 0.1 and a = 0.001 respectively. 

Both plots verify the predictions of the long-wavelength 
theory when a is small: as fc — > 0, {ai,/a)c — > 195/75 
(equation I46II . when 7 = 7/5. Also, for sufficiently small 
Q, the scales divide into the three regimes we explored ear- 
lier (Fig. 8). We find that for small a and fc ~ 1 the marginal 
curves converge to a limiting curve, differing very little from 
that appearing in Fig. 8. This confirms that the stress in this 




0.5 1 1.5 2 0.5 1 1.5 2 



Figure 6. The u'^ and u'^ vertical structure for a = 0.001, 
k = 0.1, 7 = 7/5 and a^, = X = Y = 0. The eigenfunction 
resembles the f° mode of the inviscid problem, and decays for 
these parameters. 




10"^ 10"' 10° 

k 



Figure 7. Curve of marginal viscous overstability in the (fc, a^,/a) 
plane for a = 0.1 and 7 = 7/5. The region above the curve 
is stable. As predicted by the long-wavelength theory, the curve 
approaches 197/75 ^ 2.627 as A; ^ 0. 

regime behaves like a regular perturbation upon the inviscid 
problem. 

What is interesting is the relative difficulty overstability 
encounters in the third regime, fc ~ 1, when bulk viscosity 
is added. An Qb a little larger than a is enough to quench 
overstability, whereas longer waves require a value more than 
twice a. We conclude that overstability on shorter scales 
is relatively fragile: a smaller Qb will switch it off (as will 
increasing 7). This is despite the fact that unstable modes 
in this regime will grow far more vigorously than those on 
the longer scales. The simulations of Kley et al. (1993) reflect 
this behaviour. The largest radial lengthscale resolvable in 
their model was approximately lOH, which falls comfortably 
in the intermediate- to short-scale regime of our analysis. For 
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Figure 8. Curve of marginal viscous overstability in the (fc, a^,/a) 
plane for a = 0.001 and 7 = 7/5. The region above the curve is 
stable. 

sufficiently small a, intermediate- and short-scales will be 
stable for appropriate 7 and/or a^,. But, as a increases and 
we leave the low-viscosity regime, these scales may become 
overstable, as Kley et al. report. 



6 SUMMARY AND DISCUSSION 

We have investigated the growth or decay rate of the fun- 
damental mode of even symmetry in a three-dimensional 
viscous accretion disc. Our numerical model resolves the ver- 
tical structure of the disc and its modes, treating radiative 
energy transport in the diffusion approximation. Oversta- 
bility can occur, in principle, for radial wavelengths longer 
than a few times the disc thickness. It results almost entirely 
from the fact that the unperturbed disc necessarily has a 
shear stress in order to facilitate accretion, and the varia- 
tion of this stress in the wave motion allows energy to be 
fed from the differential rotation into growing oscillations. 
This tendency competes with the viscous damping of the 
mode. Smaller contributions to the growth rate arise from 
non-adiabatic effects related to the e- and K-mechanisms in 
stars, and thermal diffusion. 

The pattern of behaviour in a low-viscosity disc (a <^ 1) 
is as follows. For radial wavenumbers kH ~ 1 (i.e. radial 
wavelengths a few times the disc thickness) the structure of 
the f° mode resembles its counterpart in an inviscid disc. It 
has a non-trivial vertical structure dictated by the stratifi- 
cation of the disc; the horizontal velocities involve vertical 
shear and are accompanied by vertical motion. The mode 
can be overstable in this short-wavelength regime, depend- 
ing on the parameters of the disc, but the presence of vertical 
shear and vertical motion acts against this tendency. 

For very long radial wavelengths {kH < a) the mode is 
dominated by horizontal motion and the vertical shear is 
removed by the action of viscous stresses between different 
strata in the disc. This type of motion is optimal for exciting 
the overstability but the (complicated) analytical criterion 
we obtain in this regime differs from the prediction of a 



two-dimensional model because the disc is not in vertical 
hydrostatic or thermal balance. 

On intermediate wavelengths (a <^ kH <g 1) a transi- 
tional behaviour is found. Here the vertical shear developing 
in the horizontal velocity perturbations leads to an enhanced 
viscous damping and the overstability is suppressed. 

In order to bring out this subtle behaviour, we have re- 
stricted our attention in this paper to viscous models of 
accretion discs based on the Navier-Stokes equation and 
the alpha viscosity prescription. We fully appreciate that 
stress in accretion discs is likely to arise from magnetohy- 
drodynamic t urbulence and to ex hibit a more complicated 
rheology (e.g. IOgilvi3l200ll . l2003l) . In particular, the non- 
zero relaxation time of the stress is known to act against 
the overstability llpg ilvic 200 J as it does in planetary rings 
iLatter fc OeilvielboOGfl . In this paper we have employed a 
bulk viscosity as a convenient surrogate for any such stabi- 
lizing tendencies. 

A major application of this work is to the theory of ec- 
centric accretion discs, since a small eccentricity can be re- 
garded as a wave of azimuthal wavenumber m — 1 and long 
radial wavelength propagating in a circular disc. Indeed, it is 
only for m = 1 waves in (nearly) Keplerian discs that a long 
radial wavelength can be maintained over an extended radial 
distance. Overstability appears as a n egative diffusi on coef- 
ficient for eccentricity in the theory of lOgilvi'3 ||200J) . which 
uses approximations similar to the long-wavelength regime 
described above. The approximations break down for inter- 
mediate wavelengths in a very low- viscosity disc, because the 
viscous stresses are inadequate to couple the different strata 
and enforce uniformity of the eccentricity with height. 

While the overstability may be suppressed by stress re- 
laxation effects, as suggested bv lOgilvig i200lf) . an alterna- 
tive possibility is that overstability may indeed occur on very 
long radial wavelengths, while shorter wavelengths require a 
treatment similar to that carried out here. 

In protoplanetary discs the stress is of uncertain origin, 
but values in the range lO"* <a< 10"^ are usually quoted 
based on observed accretion rates and disc lifetimes, while 
the angular semithickness is in the range 0.05 < iif/r < 0.1. 
In this case the long-wavelength analysis has no application, 
since the relevant wavelengths greatly exceed the size of the 
disc. Instead, the regimes of short and intermediate wave- 
lengths describe the behaviour of global eccentric modes. 
This means that a correct determination of the the damp- 
ing rate of eccentricity in a protoplanetary system requires 
a three-dimensional treatment of the disc, allowing for verti- 
cal shear in the horizontal velocities. We find that a typical 
decay rate of a global eccentric mode (with a radial wave- 
length of 20H) is 10"^ ^"^ for a = 10"^ and 2 x 10"*^ 
for a = 10"'*. These are only rought estimates as we have 
not considered the opacity regimes relevant to protoplane- 
tary discs. For a = 10~^ global modes may be overstable. 

In other accretion discs H/r may be much smaller and 
a somewhat larger. In this case the long-wavelength regime 
applies to global eccentric modes. The fact that the oversta- 
bility is suppressed for intermediate wavelengths means that 
only the largest-scale modes may be permitted to grow. This 
would neatly account for the occurrence of global eccentric 
modes in decretion discs around Be stars. 

No modes other than the mode were found to exhibit 
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overstability. In particular, a similar mechanism would not 
permit global warping of discs to be excited. 
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where starred quantities are dimensionless and satisfy the 
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dF, 
di7 
dJ\ 
dz, 

r-H, 

p, dz, = 1, (AlO) 



(A6) 

(A7) 

(A8) 
(A9) 



together with the boundary conditions p,[±Ht) = 
Tt{±Ht) = 0. For reasonable values of X and Y , this prob- 
lem has a unique solution that can be obtained numerically, 
with determined as an eigenvalue. 

Near the upper surface, the limiting behaviour is of 
the form K. constant, T, oc (_ff* — z*), p* oc — 



\(3-Y)/{l+X) 



p, OC (i?, — Zf 



The vertically averaged kinematic viscosity v is given 



by 



11 dz 



H 



aUpUn 



p, dzt, 



and therefore 

vY, oc j](io+3x-2y)/(6+x-2y)^ 



(All) 
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APPENDIX A: VERTICAL STRUCTURE OF 
THE DISC 

Dimensional analysis of the equations of vertical structure 
leads us to introduce characteristic units U for the variables 
H , p, p, T and F according to 



16a 



Uf = aAA^n-^UHUp. 



(Al) 

(A2) 
(A3) 
(A4) 
(A5) 



APPENDIX B: LONG- WAVELENGTH LIMIT 

We consider the limit fc of equations 1251 - 12911 . One 
possible scaling is 

s ^ k'^S2 + 0{k^), 



it^c = ku'^i + O(fc^), 

= ku'yl + 0{k^), 

u'^ = k^u^2 + 0{k^), 
P ^ p'o + O(fc^), 
p' =p'o + 0{k^), 
which yields 

- 2Qpu'yi = -ipo + dzipdzu'^i), 
2{Q - A)pu'^i = -2iyl/io 4- dz{p.dzUy-^), 
= p'ogz - dzp'o, 

S2p'o -f u'^2dzp + p{iu'^i + dzu'^2) = 0, 
0= (7-l)(4AVo-9zFio)- 



Therefore 
S2 / p'odz 



pu'xi dz, 



pu',1 dz = -2iA / 



fio dz. 



(Bl) 
(B2) 
(B3) 
(B4) 
(B5) 
(B6) 

(B7) 
(B8) 
(B9) 
(BIO) 
(Bll) 

(B12) 
(B13) 



Since the hydrostatic and thermal balances are maintained 
in the perturbed equations to leading order. 
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/: 



-H 

and so 

S2Sd = 



9E 



■2^0- 



n - A / 9E 
The dispersion relation for this mode is therefore 

Another possible scaling is 
I 

J-y2 - 



+ k\'2 + O(k^) 



< = ku',1 + 0(k^), 
p' = kp'^+0{k'), 
p = kp[+0{k^). 



(B14) 
(B15) 

(B16) 

(B17) 
(B18) 
(B19) 
(B20) 
(B21) 
(B22) 



The horizontal components of the equation of motion at 
leading order yield 



p [sou'yo + 2{Vl - A)u'^q\ = dz{pdzUyo). 



(B23) 
(B24) 



If M^oi depend on z, we obtain viscously damped 

epicyclic motions at this order. The modes satisfy a Sturm- 
Liouville equation and only one mode is undamped at lead- 
ing order, being independent of z. For this solution, the 
right-hand sides vanish and 

So + 4n(n - yl) = 0. (B25) 

Therefore so = ±iQ,r, u'^^ = U (independent of z) and li^o = 
V = (so/2r2)L''. The remaining equations at leading order 
yield 

+9. [2/ia.ii;i + (Mb - lp){iU + d.u',^)] , (B26) 
sop'i + u',id,p + p{iU + d,u',^) = 0, (B27) 



sopi -I- u';,idzP + 7p(i(7 + i9zM^i) 

= (7 - l)(4AVi - iA^^V - 9.Fii), 



(B28) 



Ti 
T 

I 

Ml = 



V 



Sup 

+ 

El 
p ' 



(i + x)^ + (3-y)^ 



api 



(B29) 

(B30) 
(B31) 



These can be solved, in principle, for u'^i, p'l, p'l, T[, F!,i (see 
below). Then the horizontal components of the equation of 
motion at order k yield 

p(s2(7 -f- sou^2 - 2^11^2) = -ip'i 

+i [2fiiU + (pb - iM)(iC/ + d,u',r)] 

+dz [^(iMzi + 9^14^2)] , (B32) 



P [s2V + sou'y2 + 2(n - A)u'^2] = -2iA^i'i - iiV 

+d,{pd,u'y2). (B33) 

We eliminate u'^2 and u'y2, obtaining the solvability condi- 
tion for these equations, by taking so times equation 1B321 
plus 2f2 times equation (IB33L then integrating vertically: 



j I So [-ip'i - (Mb + \p^)U 



+ (Mb — |M)i9zit'2i] — 4f2Ai/i'i > d 



(B34) 



Note the contribution to this equation from the vertical mo- 
tion, which also changes the relation between p'^ and U . 

Under convenient assumptions [R, 7, a, ctb b eing in- 
dependent of z) we have the analytical solution (cf. lOeilviel 
.200 It) 

(B35) 
(B36) 
(B37) 
(B38) 
(B39) 



u'zl = SqWz, 

p'l = -Wzdzp + App, 
p'l = - WzdzP + App, 
T[ = -Wzd.T + ArT, 
F^i = -Wzd.F, + ApF, 



In this solution, the Lagrangian perturbations of all quanti- 
ties are proportional to the unperturbed quantities, meaning 
that the disc undergoes a kind of homogeneous (but dynam- 
ical) expansion. We find 

slW = nl{~2W ~Ap + Ap) - {at, + ^a)QlW{so/ft) 

-i(ab-iQ)n'([//fi), (B40) 



soAp = -if/ - soW, 



(B41) 



soAp + 7(1(7 + soW) 

= (7 - 1) [4A''{a/n){Ap + W-Af) 

-4A{a/Q.)iV] , (B42) 

Af^-{1+ X)Ap + {4^ Y){Ap ^Ap)^W. (B43) 

We solve for W , then return to find S2 using 

H nH 

p'iAz^{W + Ap)P, P= / pAz. (B44) 

H J -H 

Therefore 

/ 4Aa\ 
2(E/P)s2(7 = -i M + — j {W + Ap) 

-(ab + |a)(t//n) + (at - |a)i(so/J7) (B45) 

From the above equations S2 can be determined alge- 
braically. 



